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Abstract: A procedure for counting the number of eigenvalues of a matrix in a 
region surrounded by a closed curve is presented. It is based on the application of 
the residual theorem. The quadrature is performed by evaluating the principal 
argument of the logarithm of a function. A strategy is proposed for selecting 
a path length that insures that the same branch of the logarithm is followed 
during the integration. Numerical tests are reported for matrices obtained from 
conventional matrix test sets. 
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Denombrement de valeurs propres 
dans le plan complexe 

Resume : Nous proposons une approche pour compter le nombre de valeurs 
propres d'une matrice, situees dans un domaine borne du plan complexe. L'approche 
est fondee sur l'application du theoreme des residus. reintegration nous ramene 
a revaluation de Fargument principal du logarithme d'une fonction. Nous 
proposons une strategie pour le choix du pas qui permette de rester sur la meme 
branche du logarithme pendant l'integration. Des resultats numeriques sont 
obtenus a partir de tests conduits sur des matrices tirees d'ensembles classiques 
de matrices. 

Mots-cles : Valeurs propres, resolvante, determinant, logarithme complexe 
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1 Introduction 

The localization of eigenvalues of a given matrix A in a domain of the complex 
plane is of interest in scientific applications. When the matrix is real symmetric 
or complex hermitian, a procedure based on computations of Sturm sequences 
allows to safely apply bisections on real intervals to localize the eigenvalues. 
The problem is much harder for non symmetric or non hermitian matrices and 
especially for non normal ones. This last case is the main concern of this work. 
Proceeding by trying to compute the eigenvalues of the matrix may not always 
be appropriate for two reasons. 

First most of the iterative methods frequently used to calculate eigenvalues 
of large and sparse matrices may loose some of them, since only a part of the 
spectrum is computed, and as such there is no guarantee to localize all the 
eigenvalues of the selected domain. When a shift-and-invert transformation 
is used, the eigenvalues are obtained in an order more or less dictated by their 
distance from the shift, and if one eigenvalue is skipped, there is no easy strategy 
that allows to recover it. 

Second the entries of the matrix may be given with some errors and then 
the eigenvalues can only be localized in domains of C. 

Many authors have defined regions in the complex plane that include the 
eigenvalues of a given matrix. One of the main tool is the Gershgorin theo- 
rem. Since a straight application of the theorem often leads to large disks, some 
authors extended the family of inequalities for obtaining smaller regions by in- 
tersections which include eigenvalues (see e.g., [7JIH])- Other techniques consist 
to consider bounds involving the singular values (see e.g., [I]), the eigenvalues 
of the hermitian part and the skew- hermitian part of the matrix (see e.g., [5]), 
or the field of values of inverses of the shifted matrices (see e.g., |10|V 

For taking into account, possible perturbations of the matrix, Godunov[S] 
and Trefethen [14) have separately defined the notion of the of e-spectrum or 
pseudospectrum of a matrix to address the problem. The problem can then be 
reformulated as that of determining level curves of the 2-norm of the resolvent 
R(z) = {zl - A)^ 1 of the matrix A. 

The previous approaches determine a priori enclosures of the eigenvalues. A 
dual approach can be considered: given some curve (r) in the complex plane, 
count the number of eigenvalues of the matrix A that are surrounded by (r). 
This problem was considered in [5] where several procedures were proposed. 
In this paper, we make some progress with respect to the work in (5j. Our 
procedure is based on the application of the residual theorem: the integration 
process leads to the evaluation of the principal argument of the logarithm of the 
function g(z) = dct((z + h)I — A)/Act(zl — A). This function is also considered 
in [S] to count the eigenvalues when a nonlinear eigenvalue problem is perturbed. 

This work is mainly concerned with the control of the integration path so 
as to stay on the same branch along an interval when evaluating the principal 
argument of a logarithm. 

In section [5J we present the mathematical tools. In section [3J we present 
the basis of our strategy for following a branch of the logarithm function and 
conditions for controlling the path length. Section 0] deals with the implemen- 
tation of our method: we show how to safely compute the determinant and how 
to include new points along the boundary. In section [S] we present numerical 
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test results carried out on some test matrices and in section we conclude with 
some few remarks and future works. 

2 Mathematical tool and previous works 

In this section we present the Cauchy's argument principle and some previous 
works on counting eigenvalues in regions of the complex field. 

2.1 Use of the argument principle 

The localization of the eigenvalues of matrix A involves the calculation of deter- 
minants. Indeed let (T) be a closed piecewise regular Jordan curve (piecewise 
C 1 and of winding number 1) in the complex plane which does not include eigen- 
values of A. The number Nr of eigenvalues surrounded by (r) can be expressed 
by the Cauchy formula (see e.g., [T21 113| ): 

where f(z) = det(zl — A) is the characteristic polynomial of A. 

If 7(i)o<t<i is a parametrization of T the equation ([T]) can be rewritten as 

^ = ^f^§-i'm. (2) 

2"r Jo A7W) 

The primitive ip defined by 

is a continuous function which is a determination of log(/ o 7) (e.g. see |13|): 

log/(7(*)) = i^|/(7(*))l+iar 5 (/(7(t))), f e [0,1]. 
It then follows that 

Z7T 

where </?/(l) is the imaginary part of tp(l) since its real part vanishes. 

2.2 Counting the eigenvalues in a region surrounded by a 
closed curve 

In [S], two procedures were proposed for counting the eigenvalues in a domain 
surrounded by a closed curve. 

The first method is based on the series expansion of log(J + hR(z)), where 
R(z) = (zl — A) -1 , combined with a path following technique. The method uses 
a predictor - corrector scheme with constant step size satisfying the constraint 

\ipi{z + Az) - tpi(z)\ < it, 
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for a discrete list of points z. The implementation of the algorithm requires the 
computation of a few of the smallest singular values and the corresponding left 
and right singular vectors of (z I — A) ; they are used to follow the tangent to 
the level curve of the smallest singular value of (zl — A) . 

In the second procedure, the domain is surrounded by a parameterized user- 
defined curve z = j(t) and thus 

1 /•T(i)|det( 7 (t)/-A) 



NT ~2i^l {0) det( 7 (t)I-A) dt (3) 
Since 7(0) = 7(1), the function j(t) defined on [0, 1], can be extended onto K 

by 

7ext(t) = j{t mod 1). 

By subdividing the interval [7(0), 7(1)] into subintervals of equal length, and by 
assuming that "f ext £ C°° , they make use of a fundamental result from quadra- 
ture of periodic function to prove an exponential convergence of the integral. 
The method is compared to other integrators with adaptive step sizes. 
Each of these methods makes use of the computation of 

detmi-A) 



\det( 7 (t)I - A)[ 



which is efficiently computed through a LU factorization of the matrix {"~f{t)I — 
A) with partial pivoting. In order to avoid underflow or overflow, the quantity 
is computed by 

det(j(t)I - A) 

\det(>y(t)I - A)\ ~ 11 K| 

where un is the i-th diagonal element of U in the LU factorization. The product 
is computed using the procedure that will be described later on in section @] 

Our work, which can be viewed as an improvement of jS], mostly deals with 
the control of the integration so as to stay on the same branch along an interval, 
during the evaluation of the principal argument of the logarithm of the function 
g(z) defined in the introduction. 



3 Integrating along a curve 

In this section, we describe strategies for the integration of the function g(z) = 
^y-j , where f(z) = det(zl — A), along the boundary of a domain limited by a 
user-defined curve (T) that does not include eigenvalues of A. 

3.1 Following a branch of log(f(z)) along the curve 

To simplify the presentation and without loss of generalization, let us assume 
that r = {J^Lq 1 [zi, Zj+i] is a polygonal curve. 

Let Arg(z) G (— 7r,7r] denote the principal determination of the argument of 
a complex number z, and arg(z) = Arg(z) (2ir), be any determination of the 
argument of z. In this section, we are concerned with the problem of following 
a branch of log(/(z)) when z runs along (T). The branch (i.e. a determination 
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arg of the argument), which is to be followed along the integrating process, is 
hxed by selecting an origin zq £ (T) and by insuring 

arg (/(zo)) = Arg(/(z )). (4) 
Let z and z + h two points of (r). Since 

(z + h)I-A = (zI-A) + hI 

= {zl - A)(I + hR(z)), 

where R(z) = (zl — A) -1 , it then follows that 

f(z + h) = f(z)det(I + hR(z)). (5) 

Let $ z (h) = det(I + hR{z)), then 

' +k tMdz = log(f(z + h))-log(f(z)) 
f .f(z + h) 

= log($ z (ft)) 

= log\$ z (h)\+i&rg(® z {h)). 

In the previous approach [3], given z, the step h is chosen such that condition 

|arg(*,(Ji))| < 7T, (6) 

is satisfied. In [5] condition ^ is only checked at point z + h but we want 
the condition to be satisfied at all the points s G [z, z + h], so as to guarantee 
that we stay on the same branch along the interval [z, z + h]. We need a more 
restrictive condition which is mathematically expressed by the following lemma: 

Lemma 3.1 (Condition (A)) Let z and h be such that [z,z + h] C (T). 

|Arg(* r («))| <tt, VsG [0,h], (7) 

then, 

arg (/(z + h)) = arg (/(z)) + Arg(* a (fc)), (8) 

where arg is the determination of the argument determined as in by an a 
priori given zq € (r). 

Proof. We prove it by contradiction. Let us assume that there exists k £ 
Z \ {0} such that 

arg Q (f(z + h)) = arg (f(z)) + Arg{<$> z {h)) + 2kn. 

By continuity of the branch, there exists s £ [0,h] such that \Arg($ z (s))\ = it, 
which contradicts Q. o 

Condition (O is called Condition (A). 
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3.2 Step size control 

In our approach, given z, the step h is chosen such that condition of Lemma 
(|3.1j) is satisfied. Condition (A) is equivalent to 

<f> z (s) $ (-oo,0],Vs e [0,h]. 

In order to find a practical criterion to insure it, we look for a more severe condi- 
tion by requiring that $z(s) G ^, where Q is an open convex set, neighborhood 
of 1, and included in il C C \ (— oo,0]. Possible options for f2 are the positive 
real half-plane, or any disk included in it and centered in 1. 
Since $ 2 (0) = 1, let 

= 1 + 5, with S = pe ie . 
A sufficient condition for ([7]) be to satisfied is p < 1, i.e. 

|*«(»)- 1| < 1, Vse [0,h] (9) 

This condition will be referred to as Condition (B), and, when only verified 
at z + h, i.e. 

|$,W-i|<i, (io) 

it will be referred to as Condition (B'). This last condition is the condition 
used in [S]. It is clear that Condition (B) implies Condition (A) whereas 
this is not the case for Condition (B'). 

Since it is very difficult to check ([9]), we apply the condition on the linear 
approximation ^ z (s) = 1 + s$ 2 (0) of $ z (s) at 0. Replacing function $ z by its 
tangent ^ z in leads to 

!*»(*)- 1| < i, Vse[o,h], (ii) 

which is equivalent to the following condition, referred as Condition (C): 

"" < ma- (12) 

Example 3.1 (First illustration) Let A = ^ ^ ^ J . It then follows that 

f(z) = z(z-l), 

*M = (! + £)(! + _^-), 

z z — 1 



z z — 1 

Let us assume that we are willing to integrate along the segment from z = 2 
to z = 1 + i. In order to see if intermediate points are needed to insure that 
the branch of the logarithm is correctly followed, we consider the previously 
introduced conditions on h = t(— 1 + i) where t £ [0, 1]. 

Condition (A): $2(^1) = 1 + tt + if * sa non positive real number if and 
only if h G [—2, —1] (J(— f + ilSL)- From that, it can easily be seen that the 
segment [0, — 1 + i] does not intersect the forbidden region. Therefore no 
intermediate points are needed. 
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Condition (B): this condition is equivalent to \h\\3 + h\ < 2. By studying the 
function (f>(t) = \h\\3 + h\ = y/2t\3 — t + it\, the parameter t must remain 
smaller than a ~ 0.566. 

Condition (B'): in this example, this condition is equivalent to the previous 
one, since the function <f>(t) is increasing with t. 

Condition (C): since $ 2 (' 1 ) = § + h> this condition limits the extent of the 
interval to \h\ < | or equivalently t < rj 0.471 . 

In the second example, we illustrate the lack of reliability of Condition 
(B'). 

Example 3.2 (Second illustration) Let A = XI n , where A S R and I n is the 

identity matrix of order n. It then follows that 

f{ Z ) = (z-xr, 
h 



*z(h) = (l 



z-X 

n 



z-X 



Let us assume that we are willing to integrate from z = A + 1 to z + h = A + e 10 . 
We consider the previously introduced conditions on h. 

Condition (A): |0| < 

Condition (B): \0\ < ^. 

Condition (B'): cos n9 > | which is satisfied for values that violate Condi- 
tion (A). 

Condition (C): ||| < arcsin^;, which is more severe than \0\ < ^ and there- 
fore guaranties Condition (A). 

In this example, if (r) is the circle with center A and radius 1, the step size 
must be reduced in such a way that more than 2n intervals are considered to 
satisfy Condition (A), or even 6n and 2irn intervals with Condition (B) and 
Condition (C) respectively. 

Practically, we consider that Condition (C) implies Condition (A), as 
long as the linear approximation is valid. Problems may occur when $ 2 vanishes. 
The following example illustrates such a situation. 

Example 3.3 (Critical situation) Let us consider the matrix of Examvle \3.1\ 

For z = 1/2, < f ) i/2(/i) = 1 — 4/i 2 ; and & 1 / 2 (®) = 0; an( ^ the conditions become 

Condition (A): h £ R or \h\ < 1/2, 
Condition (B): \h\ < 1/2, 
Condition (B'): \h\ < 1/2, 
Condition (C): is satisfied for all h € C. 
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4 Implementation 

In this section, we describe the numerical implementation of our method. Strate- 
gies for including new points and a procedure for safely computing the determi- 
nants are given. 

4.1 Avoiding overflows and underflows 

The implementation of our method requires the computation of 

_ det((z + h) - A) 
4j ~ det(zI-A) ■ 

In order to avoid underflow or overflow, we proceed as follows: 

For any non singular matrix M £ C nx ", let us consider its LU factorization 
PM = LU where P is a permutation matrix of signature a. Then det(M) = 
a n"=i( M ") where un £ C are the diagonal entries of U. If the matrix M is not 
correctly scaled, the product n"=i( u «) mav generate an overflow or an under- 
flow. To avoid this, the determinant is characterized by the triplet (p, K, n) so 
that 

det(A) = pK n (13) 



where: 



(p € C with \p\ = 1), and 



K = 



\ »=1 



The quantity K is computed through its logarithm: 

n 

log(JQ = -Vlog(M). 



By this way, the exact value of the determinant is not computed, as long as the 
scaling of the matrix is not adequate. 

In section [31 it was indicated that our algorithm will heavily be based on 
the computation of & z (h) = det(I + hR(z)). For h of moderate modulus, the 
determinant does not overflow. This can be verified since 



det((z + h)I - A) K%p 2 



det(z/ - A) 



4' Pi' 



where det(z/ — ^4) and dct((z + h)I — A) are respectively represented by the 
triplets (fii,Ki,n) and {pi,Ki,n). Before raising to power n, to protect from 
under- or overflow, the ratio K^jKx must be in the interval 



'Mfi 



fi\ 



where Mfi is the largest floating point number. When this situation is violated, 
intermediate points must be inserted between z and z + h. 
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4.2 Estimating the derivative 

An easy computation shows that the derivative $ z (0) can be expressed by : 

$1(0) = trace(i?(z)). (14) 

The evaluation of this simple expression involves many operations, as we show 
it now. By using the LU-factorization P(zl — A) = LU which is available at z, 
and by using (fT4")h we may compute & z (0) = Y^i=i where U = L~ 1 e i and 
Ui = (U*)~ 1 ei, with e, being the i-th column of the identity matrix. When A 
is a sparse matrix, the factors L and U are sparse but not the vectors Ui and 
Therefore, the whole computation involves 2n sparse triangular systems. Ex- 
periments showed that they involve more operations than the LU factorization 
of the matrix. Approximations of the trace of the inverse of a matrix have been 
investigated. They involve less operations than use of the LU factorization but 
they are only valid for symmetric or hermitian matrices [31 [2] • 

If the derivative is approximated by its first order approximation, sparsity 
helps. More specifically, given z and z + h, the derivative of Q z at is estimated 
by 

s 

where s = ah with a = min(10 _6 /i/|/i|, 1), and ji = max z6 r \z\- Therefore, the 
computation imposes an additional LU factorization for evaluating the quantity 
$ z (s). It is known that, for a sparse matrix, the sparse LU factorization involves 
much less operations that its dense counterpart. 



4.3 Test for including new points 



In this subsection, we describe a heuristic procedure for including new points in 
the interval [z, z + h]. In section[3J we introduced Condition (B) which is more 
severe than Condition (A) but might be easier to verify, and we proposed to 
test its linear approximation called Condition (C). Unfortunately Example 13. 31 
has exhibited that Condition (C) may be satisfied while Condition (B') and 
therefore Condition (B) is violated. To increase our confidence in accepting 
the point z + h, we simultaneously check Condition (C) and Condition (B'). 

When Condition (C) is violated, we insert M regularly spaced points be- 
tween z and z + h where 



M = min 



\h\ \K 



(15) 



with M max being some user defined parameter. 

In addition, we insist that Condition (C) is satisfied at each bound of the 



segment \z, z + h]. Therefore, on exit, the condition \h\ < , > , s , must also be 

11 l*,+h(0)l 

guarantied. When it is violated, we insert the point z + h/2 in the list. 
The following example illustrates the effect of this step size control. 



Example 4.1 Let A be the random matrix 



A 
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/ 


-0.63 


0.80 


0.68 


0.71 


-0.31 \ 




-0.81 


0.44 


-0.94 


0.16 


0.93 




0.75 


-0.09 


-0.91 


-0.83 


-0.70 




-0.83 


-0.92 


0.03 


-0.58 


-0.87 


V 


-0.26 


-0.93 


-0.60 


-0.92 


-0.36 ) 
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1.5r 



0.5 



-0.5 



-1.5 



-1.5 -1 -0.5 0.5 1 1.5 



Figure 1: Example 14.11 The eigenvalues are indicated by the stars. The polyg- 
onal line is defined by the 10 points with circles; the other points of the line 
are automatically introduced to insure the conditions as specified in section [ 
(M max = 1 in CE5])). 



The polygonal line (T) is determined by 10 points regularly spaced on the circle 
of center and radius 1.3. In Figure]^ are displayed the eigenvalues of A, the 
line (r) and the points that are automatically inserted by the procedure. The 
figure illustrates that, when the line gets closer to some eigenvalue, the segment 
length becomes smaller. 



4.4 Global algorithm 

The algorithm is sketched in Table[TJ From a first list Z of points, it extends the 
list Z in order to determine a safe split of the integral (fTJ) . The complexity of the 
algorithm is based on the number of computed determinants. For each z £ Z , 
the numbers det{zl — A) and $ z (0) are computed; they involve two evaluations 
of the determinant. Therefore, for N final points in Z ', the complexity can be 
expressed by: 

C = 2C LU N, 

where Clu is the number of operations involved in the complex LU factorization 
of zl - A. 

When the matrix A is real and, assuming that the polygonal line (r) is 
symmetric w.r.t. the real axis and intersects it only in two points, half of the 
computation can be saved since 




where (r + ) is the upper part of (r) when split by the real axis, and T(Z) denotes 
the imaginary part of Z . 
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Algorithm : Eigencnt 

Input 

Z={edges of (r)} ; 

Mpt s = maximum number of allowed points ; 
M max = maximum number of points to insert simultaneously ; 
Output 

neg = number of eigenvalues surrounded by (T) ; 
Status (Z)=-l ; 

while Status(Z)^ and length(Z) < M pts , 
for z G Z such that Status(z)==-1, 

Compute det(zl — A) and & z (0) ; 

Status(z) = 1 ; 
end 

for z £ Z such that Status(z)=l, 

if Condition (C) not satisfied at z, 

Generate M points Z as in (|T5)) ; 

Z=Z\JZ; Status(Z)=-l ; 
elseif Condition (B') not satisfied at z + h ; 

Z=ZU {z + h/2} ; Status(z + h/2)=-l ; 
else 

Status(z)=0 ; 
end 
end 

if no new points were inserted in Z ; 
for z e Z, 

if Condition (C) is backwardly violated ; 
Z=ZU{z~h/2}\ Status(z - h/2)=-l : 

end 
end 
end 

Integral = J2 z ez J ^ r s( < ^ ) z(0)) ; n ^9 = round(Integral/ 2tt) ; 
Table 1: Algorithm for counting the eigenvalues surrounded by (17) . 
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_gl 1 , , 1 1 1 1 1 1 

-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0.5 



Figure 2: Spectrum of the matrix of Example 15. II 

5 Numerical tests 

The tests are run on a laptop Dell (Processor Intel Core i7-2620M CPU, clock: 
2.70 GHz, RAM: 4 GB). The program eigencnt is coded in Matlab. 

In the following tests, we describe the performances of the algorithm for three 
real matrices chosen from the set Matrix Market pQ. The maximum inserted 
points in an interval is M max = 10. When (r) is symmetric w.r.t. the real 
axis, only half of the integration is performed. The storage of the matrices is 
kept sparse (except for computing the spectra of the matrices of the two first 
examples) . 

Example 5.1 (Matrix ODEP400A) This matrix is a model eigenvalue prob- 
lem coming from an ODE. 



n 


Plli 


Spectral radius 


Spectrum included in 


400 


7 


4.00 


[-4,4.38e-4]x[-0.01,0.01] 



This matrix is of small order and its spectrum is displayed in Figure [21 

The first experiment consists to focusing on the right part of the spectrum by 
defining a regular polygon of 10 vertices; the polygon is centered in the origin, 
symmetric w.r.t. with the real axis as shown in figure [3] (only its upper part 
is drawn), and of radius R = 10~ 3 . Five eigenvalues were correctly found as 
surrounded by the polygon. Some statistics are displayed in the first line of 
Tabled 

The second experiment focuses on the bifurcation between real and complex 
eigenvalues in the neighborhood of —3.5. In the box [—4, —3.4] x [— 10~ 4 i, 10 _4 «], 
89 eigenvalues are counted (see the statistics in the second line of Table 0). The 
aspect ratio of the box is large. The refining process proceeds in 16 steps to 
produce 1519 intervals from the initial four. If the integral is computed by the 
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x 10" 3 



1 














0.5 









* 
























e 




* 


— * 


e 


0.5 






* 








-1 














-1.5 


-1 


-0.5 





0.5 


1 



x10" 3 



Figure 3: Example 1 5. II First experiment on the right end of the spectrum. 



Table 2: Statistics for Example 15. II 





Nb. of eigenvalues in (r) 


nb. of intervals 


elapsed time 


Exper. 1 


5 


25 


5.6e-2 s 


Exper. 2 


89 


1519 


1.1 s 



relation © at each step (hence even if the necessary conditions for correctness 
are not satisfied), it would only have been correct at the fifth step and after ; 
this corresponds to 825 intervals. This illustrates the loss in efficiency which is 
imposed by the constraint for a safe computation. 

Example 5.2 (Matrix TOLS2000) This matrix comes from a stability anal- 
ysis of a model of an airplane in flight. 



n 


Plli 


Spectral radius 


Spectrum included in 


2000 


5.96 xl0 b 


r= 2.44 xlO 3 


[-750,0]x[-r,+r] 



Two experiments consider the right part of the spectrum. A first box 
[—20, 0] x [75i, 125i] is not symmetric w.r.t. the real axis. Therefore, the in- 
tegration is not reduced. Eight eigenvalues are decounted. The second box 
[—20, 0] x [— 500i, 500«] is symmetric w.r.t. the real axis but it includes 542 
eigenvalues. The statistics are reported in Table [3] In Figure [5J the spectrum 
and two zooms on it are displayed. 





Nb. of eigenvalues in (r) 


nb. of intervals 


elapsed time 


Exper. 1 


8 


2611 


11.0 s 


Exper. 2 


542 


15669 


57.7 s 



Table 3: Statistics for Example 15.21 
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Entire spectrum of the matrix TOLS2000 



500 " 
400 - 
300 - 
200 - 
100 - 

- 
-100 - 
-200 - 
-300 - 
-400 - 
-500 - 



Exper. 1: Box=[-20,0] x [75i, 125i] 

Figure 4: Example 15.21 Spectrum of 
regions of experiments. 



Exper. 2: Box=[-20,0] x [-500i,500z] 
matrix (up) and zooms on the two 



Example 5.3 (Matrix E40R5000) This sparse matrix comes from modeling 
2D fluid flow in a driven cavity, discretized on a 40 x 40 grid and with a Reynolds 
number is Re = 5000. 



n 


\\A\\x 


Spectral radius 


Spectrum included in 


17281 


1.21 xlO 2 


r=65.5 (*) 


[0,20.2]x[-r,+r] 



(*) Estimated by the Matlab procedure eigs. 

This example shows the reliability of the proposed procedure. Computing 
the 6 eigenvalues of largest real part with the Matlab procedure eigs (it imple- 
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ments the ARPACK code) returns the eigenvalues {8.3713 ± 64.653i, 8.8025 ± 
64.876i, 16.203, 20.166}. By increasing the number p of requested eigenvalues, 
only a few of them converged: for instance for p = 20, only the two rightmost 
were found. Increasing even further up to p = 100, 14 eigenvalues were given 
back, including the already computed two rightmost and 12 additional ones with 
real parts belonging to the interval [12.2,12.9]. Therefore, the user is inclined to 
ask for the exact situation in this region. Defining the rectangle T = T + n T_ 
where T + = (14, 14 + 2i, 12 + 2i, 12) and where T_ is the symmetric of P + w.r.t. 
the real axis, the procedure eigencnt returns 



Number of eigenvalues in (r) 


number of intervals 


elapsed time 


116 


7986 


54 h 42 inn 



Table 4: Statistics for Example 15.31 



Actually, the right number of eigenvalues was already given before the last 
refining step with 3994 intervals. Taking into account the result of the experi- 
ment, after several tries of shifts in eigs, all the 116 eigenvalues surrounded by 
(r) were obtained by requesting p = 200 eigenvalues in the neighborhood of the 
shift a = 13.5 (elapsed time: 10.2s). 

6 Conclusion 

In this paper, we have developed a reliable method for counting the eigenvalues 
in a region surrounded by a user-defined polygonal line. The main difficulty to 
tackle lies in the step control which must be used during the complex integration 
along the line. The method is reliable but it involves a high level of computation. 
This is the price to pay for the reliability. In forthcoming works, a parallel 
version of the method will be developed and implemented. The code involves 
a high potential for parallelism since most of the determinant computations 
are independent. A second level of parallelism is also investigated within the 
computation of a determinant for matrices arising in domain decompositions. 
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